
alpha = -180:5:180;
beta = -90:5:90;

alpha_r = deg2rad(alpha);
beta_r = deg2rad(beta);

[alpha_mesh, beta_mesh] = meshgrid(alpha_r, beta_r);

tip = 0;
if tip==1 
    DN = abs(cos(alpha_mesh).*1);
else 
  beta = deg2rad([0 90 45]);
  alpha = deg2rad([0 0 0]);
  R=[0.25 -0.25 0;-0.25 -0.25 0; -0.25 0.25 0; 0.25 0.25 0 ];
  k=[0 cos(beta(2)) sin(beta(2))];
  dphi=2*pi*k*R';
  B=exp(-1i*dphi);
  
end

Dn = zeros(size(alpha_mesh));

for i=1:size(alpha_mesh)
    for j=1:size(alpha_mesh)
  alpha=alpha_mesh(i,j);
 
  beta=beta_mesh(i,j);
  
  k=[cos(alpha)*sin(beta) sin(alpha)*sin(beta) cos(beta)];
 
  B1=exp(1i*2*pi*k*R');
  Dn(i,j)=abs(sum(B*(B1')))^2;
    end
end

[x, y, z] = sph2cart(beta_mesh, alpha_mesh, Dn);

figure(1);
surf( x, y, z)
xlabel('x'); ylabel('y'); zlabel('z')

minc = min( [min(min(x)) min(min(y)) min(min(z))] );
maxc = max( [max(max(x)) max(max(y)) max(max(z))] );
xlim([minc maxc]); ylim([minc maxc]); zlim([minc maxc]);

% surf( alpha_mesh, beta_mesh, R)